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Abstract 

O While ordinary differential equations (ODEs) form the conceptual framework for 

J2) modelling many cellular processes, specific situations demand stochastic models to 

qh capture the influence of noise. The most common formulation of stochastic models 

I— ' for biochemical networks is the chemical master equation (CME). While stochas- 

^ tic simulations are a practical way to realise the CME, analytical approximations 

>■ offer more insight into the influence of noise. Towards that end, the two-moment ap- 

^ proximation (2MA) is a promising addition to the established anatytical approaches 

including the chemical Langevin equation (CLE) and the related linear noise ap- 
O proximation (LNA). The 2MA approach directly tracks the mean and (co)variance 

Q\ which are coupled in general. This coupling is not obvious in CME and CLE and 

O ignored b)' LNA and conventional ODE models. We extend previous derivations of 

^ 2MA by allowing a) non-elementary reactions and b) relative concentrations. Of- 

^ ten, several elementary reactions are approximated by a single step. Furthermore, 

•j-j practical situations often require the use relative concentrations. We investigate the 

r> applicability of the 2MA approach to the well established fission 3'east cell cycle 

model. Our analytical model reproduces the clustering of cycle times observed in 
experiments. This is explained through multiple resettings of MPF, caused by the 
coupling between mean and (co)variance, near the G2/M transition. 
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1 Introduction 



At a coarse level, cellular functions are largely determined by spatio-temporal 
changes in the abundance of molecular components. At a finer level, cellular 
events are triggered by discrete and random encounters of molecules |T]. This 
suggests a deterministic modelhng approach at the coarse level (cell function) 
and a stochastic one at the finer level (gene regulation) [21 El 31 El El El 
Ellinmi]- However, stochastic modelhng is necessary when noise propagation 
from processes at the fine level changes cellular behaviour at the coarse level. 

Stochasticity is not limited to low copy numbers. The binding and dissocia- 
tion events during transcription initiation are the result of random encoun- 
ters between molecules [4j. If molecules are present in large numbers and the 
molecular events occur frequently, the randomness would cancel out (both 
within a single cell and from cell to cell) and the average cellular behaviour 
could be described by a deterministic model. However, many subceUular pro- 
cesses, including gene expression, are characterised by infrequent (rare) molec- 
ular events involving small copy numbers of molecules [H 0]. Most proteins 
in metabolic pathways and signalling networks, realising cell functions, are 
present in the range 10-1000 copies per cell [111 [131 [H]- For such moder- 
ate/large copy numbers, noise can be significant when the system dynamics 
are driven towards critical points in cellular systems which operate far from 
equihbrium [HI [TBI HI]- The significance of noise in such systems has been 
demonstrated for microtubule formation [I8], ultrasensitive modification and 
demodification reactions [I2], plasmid copy number control [19], limit cycle 
attractor [2^, noise-induced osciUations near a macroscopic Hopf bifurcation 
[2T] . and intraceUular metabolite concentrations [22] . 

Noise has a role at all levels of ceU function. Noise, when undesired, may be 
suppressed by the network (e.g. through negative feedback) for robust be- 
haviour [21 [231 [211 [23 [211 [2Z|- However, all noise may not be rejected and 
some noise may even be amplified from process to process, and ultimately 
influencing the phenotypic behaviour of the cell [HI [HI [2H1 [211 [30] • Noise may 
even be exploited by the network to generate desired variability (phenotypic 
and cell- type diversification) [2l|3ll|32l[33l|3l|- Noise from gene expression can 
induce new dynamics including amplification (stochastic focusing) [61 [35| [36], 
bistability (switching between states) and oscillations [371 EH ISll HO], that 
is both quantitatively and qualitatively different from what is predicted or 
possible deterministically. 

The most common formulation of stochastic models for biochemical networks 
is the chemical master equation (CME). While stochastic simulations [H] are 
a practical way to realise the CME, analytical approximations offer more in- 
sights into the influence of noise on cell function. Formally, the CME is a 
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continuous-time discrete-state Markov process |l2l [33], [33] • For gaining intu- 
itive insight and a quick characterisation of fluctuations in biochemical net- 
works, the CME is usuaUy approximated analytically in difl'erent ways [331 [3S], 
including the frequently used the chemical Langevin approach [iG lliTl [38l li9] . 
the linear noise approximation (LNA) [151 EOl [SD [52] and the two-moment 
approximation (2MA) |53l[Mll53]. 

Of the analytical approaches mentioned above, we here focus on the 2MA 
approach because of its representation of the coupling between the mean and 
(co)variance. The traditional Langevin approach is based on the assumption 
that the time-rate of abundance (copy number or concentration) or the flux 
of a component can be decomposed into a deterministic flux and a Langevin 
noise term, which is a Gaussian (white noise) process with zero mean and 
amplitude determined by the the dynamics of the system. This separation of 
noise from the system dynamics may be a reasonable assumption for external 
noise that arises from the interaction of the system with other systems (like the 
environment), but cannot be assumed for internal noise that arises from within 
the system [3l[Sl[IIl[I3l[Sni[SI]- As categorically discussed in [37], internal noise 
is not something that can be isolated from the system because it results from 
the discrete nature of the underlying molecular events. Any noise term in the 
model must be derived from the system dynamics and cannot be presupposed 
in an ad hoc manner. However the chemical Langevin equation (CLE) does 
not suffer from the above criticism because Gfllespie [3S] derived it from the 
CME description. The CLE aUows much faster simulations compared to the 
exact stochastic simulation algorithm (SSA) [33] and its variants. The CLE 
is a stochastic differential equation (dealing directly with random variables 
rather than moments) and has no direct way of representing the mean and 
(co)variance and the couphng between the two. That does not imply that CLE 
ignores the couphng hke the LNA which has the same mean as the solution of 
the deterministic model. 

The merits of the 2MA compared to alternative approximations have been 
discussed in [S3l [Ml [SB] • In [HH], the 2MA is developed as an approximation of 
the master equation for a generic Markov process. In [53], the 2MA framework 
is developed under the name "mass fluctuation kinetics" for biochemical net- 
works composed of elementary reactions. The authors demonstrate that the 
2MA can reveal new behaviour like stochastic focusing and bistability. An- 
other instance of the 2MA is proposed in [321 [S3] under the names "mean-fleld 
approximation" and "statistical chemical kinetics". Again, the authors assume 
elementary reactions so that the propensity function is at most quadratic in 
concentrations. The authors evaluate the accuracy of the 2MA against the 
alternatives (such as LNA) for a few toy models. The derivation of the 2-MA 
for more general systems with non-elementary reactions is one motivation for 
the present paper. 
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The 2MA approaches referred to above assume absolute concentrations (copy 
number divided by some fixed system size parameter). In systems biology, 
however, models often use relative concentrations that have arbitrary units 
[SQllHniEHES]- In general, the concentration of each component in the system 
may have been obtained by a different scaling parameter, rather than using 
a global system size. For such models, the above mentioned approaches need 
modification. This was another motivation for our derivation in this paper. 

In the present paper we develop a compact derivation of the first two-moments, 
the mean and (co)variance of the continuous-time discrete-state Markov pro- 
cess that models a biochemical reaction system by the CME. This derivation is 
an extension of previous derivations, taking into account arbitrary concentra- 
tions and non-elementary reactions. The matrix form of our derivation allows 
for an easy interpretation. Using these analytical results, we develop our 2MA 
model of the fission yeast cell cycle which has two sets of ODEs: one set for the 
mean protein concentrations and the other set for concentration (co)variances. 
Numerical simulations of our model show a considerably different behaviour. 
Especially, for the weel' cdc25A mutant (hereafter referred simply as double- 
mutant), the timings of S-phase and M-phase are visibly different from those 
obtained for a deterministic model because of the oscillatory behaviour of 
the key regulator. Since the 2MA is only an approximation, we investigate 
its vahdity by comparing the statistics computed from the 2MA model with 
experimental data. 

The rest of this paper is organised as follows. In the first section we introduce 
the basic terminology and notation. Then the system of ODEs forming the 
2MA approach is presented. Next, we introduce an application to the fission 
yeast cell cycle model [59]. We present a 2MA model of the cell cycle, followed 
by a comparison to the experimental data and conclusions. The appendices 
contain full derivations of the 2MA model, further proofs and additional tables. 



2 Stochastic modelling of biochemical systems 

Imagine a well-mixed homogeneous cellular compartment of a fixed volume V 
at thermal equilibrium that contains molecules of s different kinds (each kind 
referred to as a chemical component or species) interacting in r distinct ways 
(each way referred to as a reaction channel or step). Since these biochemical 
reactions occur by random encounters of reactant molecules, the copy num- 
ber of a particular component present in the system at time t fiuctuates. The 
state of the cellular system is described by the s x 1 random vector N{t) 
whose ith element is the copy number iVj(t) of the ith species present in the 
system at time t. Each (time- varying) element Ni{t) is a stochastic process, 
where Ni{t) = rii means that molecules of the ith species are present in the 
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system at time t. The s x 1 vector n, with elements n^, is thus a sample (or a 
value) of the stochastic process N{t). The stochastic process is characterised 
by the (time-dependent) probability distribution P{n,t), that is the proba- 
bility of N{t) = n given a fixed initial condition A^(0) = n°. The probability 
distribution itself is characterised by its moments. 

We can describe the system state at time t by the s x 1 vector X{t) whose ith 
element is the concentration Xi{t) of the ith component. The concentration 
Xj(t) is, in general, the copy number Ni(t) divided by some fixed scahng 
parameter i?j specific to that component. In other words 

Ni{t) = QiXi{t), rii = QiXi . 

Each concentration Xi{t) is a stochastic process, where Xi{t) = Xi means that 
the concentration of the component at time t is Xj. The s x 1 vector x, 
with elements Xj, is thus a sample of the stochastic process X{t). The copy 
number and concentration (vectors) are related by 

N{t) = QX{t), n = Qx, 

where i? is the diagonal matrix with i?j being its ith diagonal element. 

Commonly, all components are scaled by a single parameter, in which case Q 
is a scalar known as the system size. A common choice for the system size 
is some multiple of the volume V of the system. For molar concentrations, 
the system size chosen is i? = A^^V" where Na is the Avogadro's constant. In 
systems biology, one often uses relative concentrations Xi where i?j is some 
fixed copy number specific to component i. The simplest case of relative con- 
centrations uses a single (maximum) copy number n^ax for all components. 
Note that our approach is developed for the general case which allows for rel- 
ative concentrations instead of assuming one global system-size H as done in 

HSIEIIESIEIIES]. 

If we assume that the molecules are well mixed and are available everywhere 
for a reaction (space can be ignored), then the probability of a reaction in a 
short time interval depends almost entirely on the most recent copy numbers 
(and not its earher values). In other words, the stochastic process N{t) of copy 
numbers is Markovian in continuous-time. Since changes in the copy numbers 
require the occurrences of reactions which are discrete event phenomena, N{t) 
is referred as a jump process. The Markov property implies that each reaction 
channel j can be characterised by a reaction propensity aj{n) defined such 
that, in state n, the probability of one occurrence of reaction channel j in a 
vanishingly short time interval of length dt is aj{n)dt. 

The transition from state n to the state determined by the jth reaction will 
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be represented by the following scheme 



n 



aj{n) 



n + 5'.,- 



where S,j is the jth column of the stoichiometry matrix S whose element Sij 
denotes the change in copy number of the ith component resulting from the 
occurrence of the j'th channel. Similarly the transitions towards state n from 
the state determined by the jth reaction can be represented by 



n 



aj(n—S,j 



-> n 



where the argument of the propensity function aj is n — S,j which is the 
assumed current state. Transitions away from state n will decrease the proba- 
bility P{n,t) while those towards state n will increase it. Since this is equally 
true for each reaction channel, during a short time interval of length At, the 
change in the probability is given by 

r r 

P{n, t+At)-P{n, t)^Yl Pii^-S.j, t)aj{n-S.j)At-J2 P{n, t)aj{n)At+o{At) 

3=1 j=i 

where o{At) represents terms that vanish faster than At as the later ap- 
proaches zero. As At approaches zero in the above system of equations, we 
are led to what is known as the chemical master equation (CME): 



d 
dt 



P{n,t) = Y. 



3=11 



aj{n — S.j)P{n — S.j, t) — aj{n)P{n, t) 



(1) 



We will switch between the two alternative notations jjj(t>(t) and ^ for any 
scalar quantity 4>{t). We will prefer the later when dependence on time variable 
is imphcitly clear. 

Since there is one equation for each state n and there is potentially a large 
number of possible states, it is impractical to solve the CME. In most cases, we 
are interested in the first two-moments: component-wise copy number means 



and the covariances 

Coy {N,{t),N,{t)) 



E 



(iV,(t) - E [N,{t)]) (iVfc(t) - E [Nk{t)]) 



between copy numbers of component pairs. These covariances form the covari- 
ance matrix in which the diagonal elements are component-wise variances. 



In the present paper, we are interested in the mean concentration vector /x(t) 
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with elements 

and the concentration covariance matrix a(t) with elements 

...w = cov(a-.w.a-.w) = ^h:«^ 

Hereafter, we leave out the dependence on time to simplify the notation, but 
include it occasionally when causing confusion. 



2.1 Continuous approximations of the jump process N{t) 



While the stochastic simulation algorithm and extensions provide a way to 
generate sample paths of copy numbers for a biochemical system, the need 
for repeating many simulation runs to get an idea of the probability distribu- 
tion in terms of its moments (mean and (co)variance) become increasing time 
consuming and even impractical for larger systems. Therefore attempts have 
been made towards approximations of the CME, the most notable being the 
chemical Langevin equation (CLE) by Gillespie [IH]. He obtained that con- 
tinuous approximation for the incremental change in copy number during a 
short interval [t,t + dt] where the interval length dt satisfies two conditions: (i) 
It is small enough that the propensity does not change "appreciably" during 
the interval, and (n) is large enough that the expected number of occurrences 
E [Zj(t + dt) — Zj(t)] of each reaction channel j during the interval is much 
larger than unity. That continuous approximation takes the form of the CLE 

N^{t + dt) - Nm = E S.ja, (N^t)) dt + j: S.,^a,{N^{t))dtAr,{t) . (2) 

i=i j=i 

Here N'^(t) denotes the continuous Markov process approximating the jump 
process N{t), and the set {A/'j(t)} are statistically independent Gaussian ran- 
dom variables each with zero mean and unit variance. The probability density 
function P^{n,t) of the continuous Markov process N'^{t) obeys the (forward) 
Fokker-Planck equation (FPE) gSl iij 

1^1". t) = t (- g A + 1 ± [a,(n)P'(„. t)] . (3) 

In effect, condition (i) allows a Poissonian approximation of Zj{t + dt) — Zj{t) 
and condition (ii) allows a normal approximation of the Poissonian. The two 
conditions seem conflicting and require the existence of a domain of macro- 
scopically inflnitesimal time intervals. Although the existence of a such a do- 
main cannot be guaranteed, Gillespie argues that this can be found for most 
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practical cases. Admitting that, "it may not be easy to continually monitor 
the system to ensure that conditions (i) and (ii) [..] are satisfied." He justifies 
his argument by saying that this "will not be the first time that Nature has 
proved to be unaccommodating to our purposes." [i6] . 

Generating sample paths of (2) is orders of magnitude faster than doing the 
same for the CME because it essentially needs generation of normal random 
numbers. See |65] for numerical simulation methods of stochastic differential 
equations such as (2). However, solving the nonhnear FPE (3) for the prob- 
ability density is as difficult as the CME. Therefore, on the analytical side, 
the CLE and the associated nonlinear FPE do not provide any significant ad- 
vantage. That leads to a further simplification referred to as the linear noise 
approximation (LNA) |111[15]. The LNA is a linear approximation of the non- 
linear FPE (3) obtained by linearising the propensity function around the 
mean. The solution of the LNA is a Gaussian distribution with a mean that is 
equal to the solution of the deterministic ODE model and a covariance matrix 
that obeys a linear ODE. This is the main drawback of LNA because, for sys- 
tem containing at least one biomolecular reactions, the mean of a stochastic 
model is not equal to the solution of deterministic ODEs, as shown next. 



2.2 Mean of the stochastic model 



The mean copy number for the ith component obeys the ODE 

jEmt)]=±S^,E[a,{Nit))] (4) 

which is derived in Appendix Al. In general, the expectation on the right of 
(4) involves involves the unknown probability distribution P{n,t). In other 
words, the mean copy number depends not just on the mean itself, but also 
involves higher-order moments, and therefore (4) is, in general, not closed 
in the mean unless the reaction propensity is a linear function of which 
is the case only for zero- and first-order reactions. Take the example of a 
first-order reaction X ^ Y with n denoting the copy number of its reactant 
and k denoting the reaction coefficient. The reaction propensity a(n) = kn 
(mass action kinetics) is linear in n. From probability theory, the expectation 
becomes E{kN) = kE{N) and thus we do not need to know the probability 
distribution for solving the ODE in the mean. Only if all reactions elementary 
and are of zero or first-order, we have exact equations for the evolution of 
mean: 

which corresponds to the ODE system for the deterministic model which treats 
the copy numbers n{t) as a continuous time-varying quantity that can be 
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uniquely predicted for a given initial condition. For systems containing second 
(and higher) order reactions, a{n) is a nonlinear function and the evolution of 
the mean cannot be determined by the mean alone. Instead the mean depends 
on higher-order moments, and hence the deterministic ODE model and the 
LNA cannot be used to describe the mean in (4). 



2.3 The 2MA approach 



The present section provides only a brief outline of the 2MA approach and we 
refer to the Appendix Al for a detailed derivation. 

An exact and closed representation of mean is not possible in general, as ev- 
ident from (4). The same is true for (co) variance and higher-order moments. 
One way to solve this problem is by repeating many stochastic simulation runs 
based on CME or the CLE, and computing the desired moments from the en- 
semble runs. An alternative is to find approximations to the exact ODEs such 
as (4) for the moments. The 2MA is one such attempt which assumes closure 
to the first two-moments: the mean and (co)variance. A scheme of chemical 
reactions or a system of deterministic ODEs is the starting point. From this 
are concluded the reaction propensities aj{n) which appear as coefficients in 
the CME describing the time derivative of the probability distribution P(n, t). 
By taking the first two-moments of the CME and subsequent simplifications 
followed by appropriate scaling, two sets of ODEs for the mean concentration 
vector ii{t) and covariance matrix a{t) are derived. This is followed by Taylor 
expansions of any nonlinear functions involving the propensity vector a(n). Ig- 
noring central moments of 3rd and order higher eventually leads to the 2MA 
system: 



Itt 
da 
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where the superscript T denotes transpose of a matrix and 



1 



'i j=i 




(7) 



1 



r 



k j=i 



i',k' L 



dxi'dxk' 



The derivation of these equations is given in Appendix Al. The effective flux 
on the right in (5) is the sum of a deterministic flux f{n) and a stochastic 
fl,ux £:/(/i,cr), the latter determined by the dynamics of both the mean and 
(co)variance. This influence of the (co)variance implies that knowledge of fluc- 
tuations is important for a correct description of the mean. This also indicates 
an advantage of the stochastic framework over its deterministic counterpart: 
starting from the same assumptions and approximations, the stochastic frame- 
work aUows us to describe the influence of fluctuations on the mean. This can 
be posed as the central phenomenological argument for stochastic modelhng. 

Note that (5) is exact for systems where no reaction has an order higher than 
two because then 3rd and higher derivatives of propensity are zero. In (6), the 
drift matrix A{fi) reflects the noise dynamics for relaxation to the steady state 
and the (Taylor approximation to the 2nd order of) diffusion matrix B{n) the 
randomness (fluctuation) of the individual events. The scaling by i? confirms 
the inverse relationship between the noise, as measured by (co)variance, and 
the system size. Note the infiuence of the mean on the (co) variance in (6). 

A deterministic model treats concentrations x{t) as continuous variables that 
can be predicted entirely from the initial conditions. Hence there is no noise 
term in the deterministic model and the ODEs reduce to i; = f{x). 

Since the 2MA approach is based on the truncation of terms containing 3rd 
and higher-order moments, any conclusion from the solution of 2MA must be 
drawn with care. Ideally, the 2MA should be complemented and checked with 
a reasonable number of SSA runs. 

In [531 El], the 2MA has been apphed biochemical systems, demonstrating 
quantitative and qualitative differences between the mean of the stochas- 
tic model and the solution of the deterministic model. The examples used 
in [531 El] all assume elementary reactions (and hence propensities at most 
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quadratic) and the usual interpretation of concentration as the moles per unit 
volume. In the next section, we investigate the 2MA for complex systems 
with non-elementary and relative concentrations. The reason for our interest 
in non-elementary reactions is the frequent occurrence of rational propensities 
(reaction rates), e.g. Michaelis-Menten type and Hill type, in models in the 
system biology literature (e.g. [66|). 



3 Fission yeast cell cycle modelling 

The growth and reproduction of organisms requires a precisely controlled se- 
quence of events known as the cell cycle |67]. On a coarse scale, the cell cycle 
is composed of four phases: the rephcation of DNA (S phase), the separation 
of DNA (mitosis, M phase), and the intervening phases (gapes Gl and G2) 
which allow for preparation, regulation and control of cell division. The cen- 
tral molecular components of cell cycle control system have been identified 

ISTliH]. 

Cell cycle experiments show that cycle times (CTs) have different patterns 
for the wild type and for various mutants |69l [70]. For the wild type, the 
CTs have more or less a constant value near 150 min ensured by a size control 
mechanism: mitosis happens only when the cell has reached a critical size. The 
value 150 min has been considered in [ISl [63], [70], [71] as the CT of an average 
WT cell (also referred to as the "mass-doubling time"). The double-mutant of 
fission yeast (namely weel' cdc25A) exhibits quantised cycle times: the CTs 
get clustered into three different groups (with mean CTs of 90, 160 and 230 
min). The proposed explanation for the quantised cycle times is a weakend 
positive feedback loop (due to weel and cdc25) which means cells reset (more 
than once) back to G2 from early stages of mitosis by premature activation of 
a negative feedback loop [TOl [71] . 

Many deterministic ODE models describing the cell cycle dynamics have been 
constructed (SHI [SD [IS]- These models can explain many aspects of the 
cell cycle including the size control for both the wild type and mutants. Since 
deterministic models describe the behaviour of a non-existing 'average cell', 
neglecting the differences among cells in culture, they fail to explain curious 
behaviours such as the quantised cycle times in the double-mutant. To account 
for such curiosities in experiments, two stochastic models were constructed by 
Sveiczer: The first model [701 [ZI] introduces (external) noise into the rate 
parameter of the protein Pyp3. The second model [71] introduces noise into 
two cell and nuclear sizes after division asymmetry. Full stochastic models 
that treat all the time-varying protein concentrations as random variables are 
reported in [JHIES]- They provide a reasonable explanation for the size control 
in wild type and the quantised CTs in the double-mutant type. Both models 
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employ the Langevin approach and hence require many simulation runs to 
provide an ensemble for computing the mean and (co)variance. However, the 
simulation results of stochastic models in |3Sl ESI EOl EH El] represent one 
trajectory (for a large number of successive cycles) of the many possible in 
the ensemble from which the CT statistics (time averages) are computed. We 
will see that the time-averages computed from the 2MA simulation are for the 
ensemble of all trajectories. 



3.1 The deterministic model 



We base our 2MA model on the deterministic ODE model for the fission yeast 
cell cycle, developed by Tyson-Novak in |59]. In this context, the cell cy- 
cle control mechanism centres around the M-phase promoting factor (MPF), 
the active form of the heterodimer Cdcl3/Cdc2, and its antagonistic interac- 
tions with enemies (Ste9,Slpl,Ruml) and the positive feedback with its friend 
Cdc25. These interactions, among many others, define a sequence of check 
points to control the timing of cell cycle phases. The result is MPF activity 
oscillation between low (Gl-phase), intermediate (S- and G2-phases) and high 
(M-phase) levels that is required for the correct sequence of cell cycle events. 
For simplicity, it is assumed that the cell divides functionally when MPF drops 
from 0.1. 

Table 1 lists the proteins whose concentrations Xj, together with MPF con- 
centration, are treated as dynamic variables that evolve according to 



Here f^{x) is the production fiux and ff{x) is the elimination fiux of ith 
protein. Note that the summands in the fiuxes f^{x) and ff{x) are rates 
of reactions, most of which, are non-elementary (summarizing many elemen- 
tary reactions into a single step). Quite a few of these reaction rates have 
rational expressions which requires the extended 2MA approach developed in 
this paper. The MPF concentration x^api can be obtained from the algebraic 
relation 

_ (Xi — X2) {Xi — Xtrim) /qx 
Xi 
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Table 1 

Proteins and fluxes. Here x denotes the vector of concentrations xi to xq. 



Index 

i 


Protein 


Production flux 

/i (^) 


Elimination flux 

fi (^) 


1 


CdclST 


kiM 


(fc^ + k'^X'i + A;^"x5) xi 


2 


preMPF 


(xi — X2) k^fee 


(A;25 + ^2 + /C^'xs + /C^"X5) X2 


3 


Ste9 


( ki+k'^x-i](l-x-i) 


J4+X3 


4 


SlplT 


•^4 +^mpf 


kQXi 


5 


Slpl 


(x4-a;5)x6 
'i' J7+a;4-a;5 


k^x^ + ksj^%^ 


6 


lEP 


J, {l~^6)^mpf 




7 


Rum It 




(/Cl2 + A;i2X8 + k'^Xrapi) X-j 


8 


SK 




ki4Xs 



where 



dM 
~dt 

•^trim 
k 

hb 
E 

G{a, b, c, d) 



pM 



2x1X7 



^wee ~^ (^wee ^wee) ^ (^^wee; ^wee-^mpf; <-^awee; Jiwee) 
^25 + (^25 ~ ^25) ^ (K253;mpf , ^25; <^a25? -^i25) 
Xi + X7 + i^diss, 

2a(i 



(10) 



bc + ad 



be + a(i)2 — 4(6 — a)ad 



Note that the cellular mass M is assumed to grow exponentially with a rate 
p, and the concentrations (xtrim, a;tf, k^ee, ^25) are assumed to be in a pseudo- 
steady-state to simplify the model. Note that we use a shghtly different nota- 
tion: p for mass growth rate (instead of p), Xtrim for Trimmer concentration 
and Xti for TF concentration. We have to emphasise that the concentrations 
used in this model are relative and dimensionless. When one concentration is 
divided by another, the proportion is the same as a proportion of two copy 
numbers. Hence, such a concentration should not be interpreted as a copy 
number per unit volume (as misinterpreted in |63]). The parameters used in 
the Tyson-Novak model [59| are listed in Table A3.1 in Appendix A3. 

The deterministic ODE model describes the behaviour of an 'average cell', 
neglecting the differences among cells in culture. Specifically, it fails to explain 
the experimentally observed clusters of the CT-vs-BM plot and the tri-modal 
distribution of CT PI EOl ED EH- 
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3.2 Feasibility of Gillespie simulations 

Ideally, we should repeat many runs of Gillespie's SSA and compute our de- 
sired moments from the ensemble of those runs. At present, there are two 
problems which this. The first problem is the requirement of elementary reac- 
tions for SSA. The elementary reactions underlying the deterministic model 
[59j are not known. Many elementary steps have been simplified to obtain that 
model. Trying to perform SSA on non-elementary reactions will lose the dis- 
crete event character of SSA. The second problem arises from the fact that the 
SSA requires copy numbers which in turn requires knowledge of measured con- 
centrations. All protein concentrations in the model are expressed in arbitrary 
units (a.u.) because the actual concentrations of most regulatory proteins in 
the cell are not known [22] • Tyson and Sveiczeip] define relative concentration 
Xi of the ith protein as Xi = rii/ i?j where i?, = CiN^V. Here Ci is an unknown 
characteristic concentration of the ith component. The idea is to make the 
relative concentrations Xi free of scale of the actual (molar) concentrations 
Ui/NAV. Although one would like to vary Cj, this is computationally inten- 
sive. This problem is not so serious for the continuous approximations such as 
CLE, LNA and the 2MA which are all ODEs and can be numerically solved. 

3.3 The stochastic model using Langevin's approach 

In [63j a stochastic model is proposed that replaces the ODE model (8) with 
a set of chemical Langevin equations (CLEs) 

which uses the Langevin noise terms: White noises Tf and T~ scaled by 

\Jft{^) ctnd \J fi {x) to represent the internal noise. The system parameter Q 
has been described as the volume by the author. As we discussed before, the 
concentrations are relative levels with different system size parameters. That 
means that concentrations are not the same as copy numbers per unit volume. 

Another stochastic model employing the Langevin's approach is reported in 
|48j which approximates the squared noise amplitudes by hnear functions: 

j^x.it) = fi {x{t)) + ^2DiXi{t)ri{t), 

where Di is a constant. The reason why the model dynamics /(x) are missing 
in this model is that the author wanted to represent both the internal and 

^ Personal communication. 



-x.(t) = /+(x(t))-/r(x(t)) + - 
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external noise by the second term on the right. 



3.4 The 2MA model 



For the cell cycle model, the flux / and the difl^usion matrix B, defined in (7), 
have elements 

MX) = /i W - fi [x), Bikix) iUy^k 

The ofl^-diagonal elements of B are zero because each reaction changes only 
one component, so that SijSkj = for i ^ k. Once these quantities are known, 
it foUows from (5) and (6) that the following set of ODEs: 

(j) (11) 

2Y,Aii{ii)aH + [Biiiii) + eB,Al^,a)] (12) 
I 

J2 [Ai{l^)(^ik + (^iiAkiil^)] i k (13) 
I 

approximates (correctly to the 2nd order moments) the evolution of component- 
wise concentration mean and covariance. See See Tables A3.2-A3.4 in Ap- 
pendix A3 for the respective expressions of the drift matrix A, the stochastic 
flux £/ and the correction-term sb added to the diffusion matrix B in (12). 

Having at hand the moments involving the eight dynamic variables x\ to xs, 
the mean MPF concentration can be shown to be approximately (correct to 
2nd order moments): 



dt 

dojA 

~df 
d(Tik 



A*mpf — l^l — 1^2 — 3;trim + 



•^trim 



1 + ^ /^2 



(14) 



for the mean MPF concentration with the understanding that Xtrim is in pseudo 
steady state (See Appendix A2 for the derivation). This expression for the 
average MPF activity demonstrates the influence of (co)variance on the mean 
as emphasised here. We see the dependence of mean MPF concentration /Xmpf 
on the variance an and covariance ai2 in addition to the means 1^1,^2 and 

•^trim ■ 
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100 200 300 400 500 600 100 200 300 400 500 600 

Time (min) Time (min) 



(a) (b) 

Figure 1. The time-courses of mass and MPF activity: (a) for the wild type, (b) 
for the double-mutant type. The 2MA predicted mean trajectories are plotted as 
solid lines and the corresponding deterministic trajectories as dashed lines. The 
difference between the two predictions is negligible for the wild type, but significant 
for double-mutant type. 

3. 5 Simulations of the 2MA model 



The system of ODEs (11)-(13) was solved numerically by the MATLAB solver 
odel5s |75]. The solution was then combined with algebraic relations (14). 
For parameter values, see Table A3.1. Since information about the individual 
scaling parameters i7j used in the definition of concentrations is not available, 
we have used Qi = 5000 for all i. This value has also been used in [63], although 
there is no clear justification. Note, however, that the 2MA approach developed 
here will work for any combination oi{Qi\. The time-courses of mass and MPF 
activity are plotted in Figure la for the wild type and in Figure lb for the 
double-mutant type. For the wild type, the 2MA predicted mean trajectories 
do not differ considerably from the corresponding deterministic trajectories. 
Both plots show a more or less constant CT near 150 min. Thus internal noise 
does not seem to have a major influence for the wild type. 

For the double-mutant type, the difference between the 2MA and determinis- 
tic predictions is significant. The deterministic model (8) predicts alternating 
short cycles and long cycles because cells born at the larger size have shorter 
cycle, and smaller newborns have longer cycles [59]. This strict alternation due 
to size control is not observed in experiments: cells of same mass may have 
short or long cycles (excluding very large cells that have always the shortest 
CT) [691 E]- This lack of size control is reproduced by the 2MA simulations: 
the multiple resettings of MPF to G2, induced by the internal noise, result in 
longer CTs (thus accounting for the 230-min cycles observed experimentally). 
Such MPF resettings have been proposed in [701 |7T] to explain quantised CTs. 
No such resetting is demonstrated by the deterministic model. 
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Figure 2. Variance an (of CcIcISt) and covariance au (between CcIcISt and 
preMPF): (a) for the wild type, (b) for double-mutant type. 

Note that the mean fi{t) of the 2MA describes the average of an ensemble of 
cells. Yet the MPF resettings observed in Figure (lb), near G2/M transition, 
introduce the required variability that explains the clustering of the cycle 
time observed in experiments. This is in contrast to the alternative stochastic 
approaches in [HI |63l [70], [71], [71] that use one sample trajectory rather than 
the ensemble average. 

How do we explain this significant effect of noise for the double-mutant on 
one hand and its negligible effect for the wild type on the other hand? If we 
look at expression (14), we see the influence of the variance au (of CdclSx) 
and covariance ai2 (between CdclSx and preMPF) on the mean MPF con- 
centration /impf- The two (co)variances are plotted in Figure 2a for the wild 
type and in Figure 2b for the double-mutant type. It is clear that the two 
(co)variances have very small peaks for the wild type compared to the large 
peaks for the double-mutant type. Note that the larger peaks in Figure 2b 
are located at the same time points where the MPF activity exhibits osciUa- 
tions and hence multiple resettings to G2. This suggest that the osciUatory 
behaviour of MPF near the G2/M transition is due to the influence of the 
oscfllatory (co)variances. This coupling between the mean and (co)variance is 
not captured by the deterministic model. 



It has to be reahsed that the above proposition requires vahdation since 
the 2MA approach ignores 3rd and higher-order moments. We cannot know 
whether that truncation is responsible for the osciUations in Figures 1 and 2, 
unless compared with a few sample trajectories simulated by the SSA. How- 
ever, as discussed before, the SSA cannot be performed (at present) for the 
model in consideration. Therefore we need to compare the 2MA predictions for 
the double-mutant type cells with experimental data. Towards that end, values 
of cycle time (CT), birth mass (BM) and division mass (DM) were computed 
for 465 successive cycles of double-mutant cells. Figure 3 shows the CT-vs-BM 



17 



Table 2 

Statistics over 465 successive cell cycles of the double-mutant type cells, predicted 
by the 2MA model, compared with experimental data, see |M1 Table 1]. 



Case 


MCT 


(TCT 


CVcT 




<7DM 


CVdm 




fBM 


(1) 


131 


47 


0.358 


2.22 


0.45 


0.203 


1.21 


0.24 


(2) 


138.8 


12.4 


0.09 


3.18 


0.101 


0.0319 


1.59 


0.0575 


(3) 


138.8 


17.6 


0.127 


3.25 


0.178 


0.055 


1.623 


0.0934 


(4) 


138.8 


23.9 


0.172 


3.32 


0.231 


0.0697 


1.657 


0.12 



(1) experimental data, (2) f2 = 5000, (3) (2 = 5200, (4) 12 = 5300. 



plot and the CT distribution for three different values {5000, 5200, 5300} of 
system size i?. 

To make this figure comparable with experimental data from [691 we as- 
sume that 1 unit of mass corresponds to 8.2 /xm cell length [71] • We can see 
the missing size control (CT clusters), in qualitative agreement with experi- 
mentally observed ones (see [Ml Figure 6] and [TQl Figure 5] for a comparison). 
There are more than four clusters, which may have arisen from the truncated 
higher-order moments. The extreme value of CT higher than 230 min suggests 
more than two MPF resettings. Furthermore, more than three modes in the 
CT distribution may have arisen from the truncated higher-order moments. 
Table 2 compares the statistics for the double-mutant type cells, computed 
with the 2MA approach, with data from [Ml Table 1]. Column 2-4 tabulate, 
for CT, the mean (iqt, the standard deviation ctct and the coefficient of varia- 
tion CVcT, respectively. The other columns tabulate similar quantities for the 
division mass (DM) and the birth mass (BM). We see that only the mean CT 
is in agreement with the experimental data. The mean values for both BM 
and DM are larger than the corresponding experimental values. The other 
statistics are much smaller the corresponding experimental values. This and 
the above plots suggest that the 2MA should be used with caution. However, 
another aspect of the cell cycle model deserves attention here. The way the rel- 
ative protein concentrations have been defined implies unknown values of the 
scahng parameters {i?,}. Since i?j = CiNaV, knowing the volume V does not 
solve the problem: the characteristic concentrations {Cj} are still unknown. 
Our simulations have chosen typical values H = {5000,5200,5300}. The cor- 
responding three pairs of plots in Figure 3 and rows in Table 2 demonstrate 
a dependence of the results on a suitable system size. There is no way to 
confirm these values. The scahng parameters could be regulated in a wider 
range in order to imporve the accuracy of our simulation, motivating future 
work for us. The conclusion is that the quantitative disagreement of the 2MA 
predictions can be attributed to two factors: 1) the truncated higher-order 
moments during the derivation of the 2MA, and (2) the unknown values of 
scahng parameters. 
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Figure 3. Cycle time behaviour over 465 successive cycles of the double-mutant cells, 
predicted by the 2MA model. (a,c,e): CT vs BM, (b,d,f): CT distribution, (a,b): 
J? = 5000, (c,d): i? = 5200, (e,f): i? = 5300. The plots are in qualitative agreement 
to experiments, see [SSI Figure 6] and [TDJ Figure 5] for a comparison. 
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4 Conclusions 



The recently developed two-moment approximation (2MA) [531 El] is a promis- 
ing approach because it accounts for the couphng between the means and 
(co)variances. We have extended the derivation of the 2MA to biochemical 
networks and established two advances to previous efforts: a) relative con- 
centrations and b) non-elementary reactions. Both aspects are important in 
systems biology where one is often forced to aggregate elementary reactions 
into single step reactions. In these situations one cannot assume knowledge of 
elementary reactions to formulate a stochastic model. Previous derivations as- 
sumed elementary reactions and absolute concentrations. However, numerous 
existing models in systems biology use relative concentrations. 

We investigated the applicability of the 2MA approach to the well established 
fission yeast cell cycle model. The simulations of the 2MA model show os- 
cillatory behaviour near the G2/M transition, which is significantly different 
from the simulations of deterministic ODE model. One notable aspect of our 
analytical model is that, although it describes the average of an ensemble, it 
reproduces enough variability among cycles to reproduce the curious quantised 
cycle times observed in experiments on double mutants. 



Acknowledgements 

In the process of preparing this manuscript Hanspeter Herzel (Humboldt Uni- 
versity), Akos Sveiczer (Budapest University of Technology and Economics) 
and Kevin Burrage (University of Queensland, Brisbane) helped in clarifying 
questions related to the cell cycle and stochastic modelling. We very appreciate 
their willingness to have discussed these issues with thus. M.U. has been sup- 
ported by the Deutsche Forschungsgemeinschaft (DFG) through grant (WO 
991/3-1). O.W. acknowledges support by the Helmholtz Alliance on Systems 
Biology and the Stellenbosch Institute for Advanced Study (STIAS). 



Appendices 

Al Derivation of the 2MA equations 

The progress of a particular reaction can be described by a quantity known 
as the degree of advancement (DA). We will write Zj{t) for the DA of the jth 
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reaction, where Zj{t) = zj means that the jth reaction has occurred Zj times 
during the interval [O.t). In the same interval the jth reaction will contribute 
a change of ZjSij molecules to the overall change in the copy number Ni of 
the ith component. Summing up contributions from all the reactions, the copy 
number can be expressed as 

Ni{t) = Ni{0) + j2SijZjit) . (Al.l) 

3=1 

Based on the definition of reaction propensity, the number of occurrences 
Zj(t + At) — Zj{t) during a short interval + Ai] has the probability distri- 
bution 

Pr[Z,{t + At)-Z,it) = z,\N{t)=n] 

'aj{n)At + o{At) if zj = 1 

1 - aj{n)At + o{At) if zj = (Al.2) 
_o(At) if Zj > 1 

where o{At) represents a quantity that vanishes faster than At as the later 
approaches zero. In effect, (Al.2) gives the conditional probability distribution, 
in state n, of the random progress (DA increment) Zj{t + At) — Zj{t) of the 
jth reaction during the time interval + At). The expected value of this 
short-time DA increment can be obtained from (Al.2) as 



E [Zj{t + At) - Zj{t) I N{t) = n] 

Z,it + At)-Z,it) = z,\Nit)=n] 

(A1.3) 



= E [^At + At) - Z,{t) = z, I N{t) = n] 

3=0 



= aj{n)At + o{At) 

which is conditioned on N{t) = n. The unconditional expectation of the DA 
increment can be obtained by summing the probabilities P{n, t) weighted by 
the above conditional expectation over all possible states n: 

E [Z,{t + At) - ZM = ^ E [Z,{t + At) - Z,{t) I Nit) = n] P{n, t) 

n 

= 5^a,(n)p„(i)Ai + o(Ai) 

n 

= E [aj(^N{t))] At + o{At) 
which for vanishingly small At leads to the ODE 

|E[Z,(t)]=EK(Ar(t))] (A1.4) 
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Thus the mean propensity of a particular reaction can be interpreted as the 
average number of occurrences (DA) per unit time of that reaction. Take the 
expectation on both side of the conservation (Al.l) to obtain 



d 
dt 



which proves (4) in the main text. It is interesting to note that the above 
ODE is a direct consequence of mass conservation (Al.l) and definition of 
propensity because we have not referred to the CME (which is the usual 
procedure) during our derivation. 



Dividing (4) by f2i gives the ODE for the component mean concentration 

d 



dt 



/i,(t) = E[/,(x(t 



(A1.5) 



where 



1 



is the total flux of component i in state x. 

Suppose the propensity %(n) is a smooth function and that central moments 
E [{N — //)"*] of order higher than m = 2 can be ignored. In that case, the 
Taylor series expansion of flux fi{x) around the mean is 



fi{x) = fiif^) + 



df^ 
dx'^ 



X=fl 



dxdx'^ 



{x-n) + 



X=fl 



Expectation of the 2nd term on the right is zero. Expectation of the 3rd term 
can be written as 



1 



dxkdxi 



x=n 



Note that the Taylor expansion in powers of x — is more convincing than 
that in powers of n — E(n) because higher-order terms vanish quicker in the 
former. Having arrived at this point, ignoring terms (moments) higher than 
2nd order, we can write: 



djij 
dt 



(A1.6) 



for mean component concentration and 

d/i 



dt 
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for the mean concentration vector. This last equation proves (5) in the main 
text. Here the term £/•(//, (t) is the internal noise that arises from the discrete 
and random nature of chemical reactions. Note that this term has been derived 
from the CME instead of being assumed hke external noise. This shows that 
knowledge of fluctuations (even if small) is important for a correct description 
of the mean. This also indicates an advantage of the stochastic framework 
over it deterministic counterpart: starting from the same assumptions and 
approximations, the stochastic framework allows us to see the influence of 
fluctuation on the mean. Note that the above equation is exact for systems 
where no reaction has an order higher than two because then 3rd and higher 
derivatives of propensity are zero. 



Before we can see how the covariance a evolves in time, let us multiply the 
CME with riiUk and sum over afl n, 

dP(n,t) ^ 
^riiUk — ^jp- = ^nink^[aj{n - S.j)P{n - S.j,t) - aj{n)P{n,t)] 



n j=l 
r 



n j=l 



{ui + Sij) {uk + Skj) aj{n)P{n, t) - ninkaj{n)P{n, t) 



= E E ('^kSij + riiSkj + SijSkj) aj{n)P{n, t) 

n j=l 

where dependence on time is implicit for all variables except n and s. Dividing 
by and recognising sums of probabflities as expectations, 

^^d^ ^ E [X,MX)] + E [XMX)] + -^TfE [B^>c{X)] 



where B{x) is the diffusion matrix with elements 

1 

Bik{.x) = —=.Y,SijSkjaj{Qx) . 
\/ iiiiik j=\ 

The relation Uik — E [XjXjt] — can be utilised to yield 



dui 



ik 



dt 



E [{Xk - i^k) MX)] + E [{Xi - i^i) fk{X)] + 



E[Bik{X)] (A1.7) 



for the covariances between concentrations of component pairs. The argument 
of the flrst expectation in (Al.7) has Taylor expansion 



fi{x) {Xk - Hk) = MlJ') {Xk - IJ'k) + 



dx'^ 



{x - n) {xk - Hk) + 



Expectation of the flrst term on the right is zero. Ignoring 3rd and higher-order 
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moments, the first expectation in (Al.7) is then 

I 

where A{x) is the drift matrix (the Jacobian of f{x)) with elements 

dU{x) 



Aikix) 



dxk 



By a similar procedure, the second expectation (Al.7) is 

I 

correct to 2nd-order moments. The element Bik{x) of the diffusion matrix has 
Taylor expansion 



dB 



ik 



dx^ 



x=ii 



d'B, 



ik 



dxdx^ 



(x -/i) + ■ ■ ■ 



x=ii 



Taking term-wise expectation, and ignoring 3rd and higher-order moments, 

E[BikiX)] = BM+eBjfi,cr) 



where 



1 



i'k' 



d^B 



ik 



dxi'dxk' 

Having these results at hand, we can now write 



Ci'k' ■ 



X = fl 



dt ^ ^ '\J~fljJJk 
for the component-wise covariances. In matrix notation 

proves (6) in the main text. The drift matrix A(/i) reflects the dynamics for 
relaxation (dissipation) to the steady state and the diffusion matrix B{fi) 
the randomness (fluctuation) of the individual events [1]. These terms are 
borrowed from the fluctuation- dissipation theorem (FDT) [76],[77j, which has 
the same form as (6). Remember that (6) is exact for systems that contain only 
zero and flrst-order reactions because in that case the propensity is already 
linear. 
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A2 Mean MPF concentration 



To find the mean MPF concentration, we start with the MPF concentration 



2^mpf = {xi -X2)(l 



•^trim \ 



-) 



X2 

I 

Xl 



The ratio ^^/xi can be expanded around the mean, 



X2 



X2 



Xl A<i 1 + 



Ml 



(Xi - /Xl) X2 , (Xi - III) X2 , 

X2 1 o h 



;/2 



Taking expectation on both sides. 



E 



X2 
Xl 



-E 



= 1e 

_ 1 



1 + 



Ml 



■ {Xi-iii)X2 , {Xi-i^iyx2 , 

A2 1 5 h 

1^2 h 2- 



Finally, the mean MPF concentration follows from the expectation of x-^apf to 
be 



A*mpf — 1^1 — 1^2 — Xxrim + 

thus proving (14) in the main text. 



A3 Parameters and coefficients of the 2MA equations 
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Table A3.1 

Parameter values for the Tyson-Novak cell cycle model of the fission yeast (wild 
type) [SH]. All constants have units min~^, except the Js, which are dimension- 
less Michaelis constants, and i^dissj which is a dimensionless equilibrium constant 
for trimer dissociation. For the double-mutant type, one makes the following three 
changes: k'^^^ = 0.3, = = 0.02 . 



ki5 = 0.03, k2 = 0.03, k2 = l, k'^' = 0.1, /cg = 1, fcg = 10, J3 = 0.01, 

k'^ = 2,k4 = 35, J4 = 0.01, k'^ = 0.005, k'^ = 0.3, fee = 0.1, J5 = 0.3, 

k'j = 1, /cg = 0.25, J7 = Js = 0.001, Js = 0.001, kg = 0.1, kio = 0.04, 

Jg = 0.01, Jio = 0.01, kn = 0.1,ki2 = 0.01, k[2 = 1, ^12 = 3, i^diss = 0.001, 

ki3 = 0.1, ku = 0.1, A:i5 = 1.5, k[e = 1, fc'/g = 2, J15 = 0.01, Jie = 0.01, 

Kwee = 0.25, Mwee = 1, Jawee = 0.01, Jiwee = 0.01, Va25 = 1; ^25 = 0.25, 

Ja25 = 0.01, Ji25 = 0.01, Ce = 0-15, Ce = 1-3, ^25 = 0-05, A^^'g = 5, p = 0.005. 



Table A3. 2 

Rows of the drift matrix A of the 2MA cell cycle model. We here use to denote 
the ith row of 8 x 8 identity matrix. 



Index i 




1 


— {k'2 + k2X3 + kl^'x^) ei — ^2^163 — /cg'^i^s 


2 


A^wccei - (AJwcc + k25 + k'2 + k2X3 + k2X5) 62 - A:2^2e3 - fc2"x2e5 


3 


(fc^a;g-|-fc4Xmpf) J4 (fc^-l-fcjj'xs) J3 {1-X3)k'^ k'^xz 

' {J3+l-xs,f ' Js+l-x^^^ Ja+x^^^ 


4 


-^664 


5 


kjJjXQ 


U_ 1 kfJ-jXQ 1 fcgJg 


1 (2^4-a;5)fc7 


(Jj + X4,-Xc,)^ 






6 


r kQX^ptJQ fcloJlO 1 

_(jg+i_^g)2 1 (Jlo-hx6)^ 


7 


- (/i;i2 + /c'i2X8 + fca'a^mpf ) 67 - k[2X7(is 


8 


-ki^es 
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Table A3. 3 



Stochastic flux, the correction-term added to the deterministic flux in (5). 



Index i 




1 




2 




3 






iJ4+x3f (Js+l-xaf 


^3 (J3+l-a:3)' {J4+X3)'' 


4 





5 


fe7J7a;6(2iT45— (T44— (T55) . kYJr(<T46—<T5e) 1 ksJa _ 

{J7+X4-Xr,)'-^ iJT+X4-Xsf {Js+X5? 


6 


r fcioJio *:9a;mpf'^9 1 
.(JlO+Xef (J9+l-X6)^ 


7 




8 






Table A3.4 



Correction-term added to Bii{x) in (12). 



Index i 


ei3«(x,c7) = lY.k,i dxkdli'^ki 


1 


k'^ai3 + k'^'ai5 


2 


^2 '^23 + fe2''''25 


3 


(fc4X8-|-fc4a;mpf)-'4 (fcg+fcgXs) J3 


k'^J:iCr35 k'^J4(TZ» 


(J4-ha:3)^ (J3-hl-X3)^ 


'^^^ {Js + l-Xif ' (J4-hX3)^ 


4 





5 


A;7J7X6(2ct45— 0-44 — 0-55) . fc7 J7((T46 — crse) fca^S - - 
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